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Abstract 

We re-explore the kinetics of spinodal decomposition in off-critical polymer 
blends through numerical simulations of the Cahn-Hilliard equation with the 
Flory-Huggins-De Gennes free energy functional. Even in the absence of ther- 
mal noise, the solution of the discretized equation of motion shows coarsening 
in the late stages of spinodal decomposition without evidence of pinning, re- 
gardless of the relative concentration of the blend components. This suggests 
this free energy functional is not sufficient to describe the physics responsible 
for pinning in real blends. 



*Permanent address: Dipartimento di Scienze Fisiche, Universita degli studi di Napoli, Pad. 19, 
Mostra d'Oltremare, Napoli 80125, Italy. Internet address: castellano@na.infn.it 
^Internet address: glotzer@ctcms.nist.gov 

1 



I. INTRODUCTION 



Experiments on spinodal decomposition in polymer blends show that the coarsening 
process may sometimes slow dramatically or even cease before reaching equilibrium |l| - 
||. In these systems, spinodal decomposition — which is the process by which a ther- 
modynamically unstable mixture demixes to a stable, phase-separated equilibrium state 
" EH ~ proceeds normally for some time following a quench to the unstable region 
(T < T s ), and then stops. A break-up of the characteristic, interconnected pattern is ob- 
served to precede this pinning phenomenon. The nonequilibrium, microphase-separated 
blend has been observed to remain in this pinned state over an appreciable time scale 
where little domain growth occurs. The eventual breakup of the evolving morphology 
into separated droplets is a natural consequence of the asymmetric composition in an 
off-critical blend ||14|| ; nevertheless, it may also occur in near-critical blends due to other 
forces (e.g. gravity). Polymer blends in which pinning has recently been observed for off- 
critical composition include X-7G/poly (ethylene teraphthalate) (a liquid crystalline poly- 
mer /ho mopolymer blend) [|IJ, poly(styrene-ran-butadiene)/polybutadiene [|J, poly(styrene- 
ran-butadiene)/polyisoprene ||], and polybutadiene/polyisoprene ||. 

The specific mechanism responsible for pinning in these blends is poorly understood, 
and is currently a topic of considerable discussion. While there is general agreement that 
growth stops soon after the breakup into separated "droplets" or "clusters" (a so-called 
"percolation-to-cluster transition" @JT^]), the mechanism that prevents further coarsening 
of disconnected domains remains to be clarified. One intriguing scenario points to an entropic 
barrier as the reason for the observed arrested growth of off-critical phase separating blends. 
Kotnis and Muthukumar (KM) Jl5| have suggested that due to the connectivity of the chains 



and the reduced conformational entropy near domain interfaces [16], the usual evaporation- 



condensation mechanism of coarsening |T7j observed in small-molecule mixtures is suppressed 
in polymer blends, and instead coarsening occurs via parallel transport of chains along the 
interface JjlJ. Consequently, KM postulate that if the concentration of the minority-rich 
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phase becomes smaller than the percolation threshold, the parallel coarsening mechanism 
will be inhibited and the clusters will "freeze" after an initial growth period. 

Hashimoto and coworkers have instead postulated that the enthalpy of mixing, rather 
than the entropy, provides the barrier to further coarsening following the percolation-to- 
cluster transition ||. They argue that the increase in enthalpy of mixing suffered upon 
removing a chain of species A and degree of polymerization N from the surface of an A-rich 
domain is AH m i x oc xNk B T, where \ is the Flory interaction parameter, k B is Boltzmann's 
constant, and T is temperature. In the strong segregation limit \N ^> 1, and thus evapora- 
tion of the chain from the domain surface, which would occur with a Boltzmann probability 
proportional to exp(— AH mix /k B T), is highly unfavorable. Thus, when the parallel transport 
mechanism is eliminated by the breakup into droplets, coarsening ceases. 

In this paper, we re-explore the kinetics of spinodal decomposition in off-critical polymer 
blends described by the Flory-Huggins-De Gennes (FHDG) free energy functional, through 
numerical simulations of the Cahn-Hilliard (CH) equation. In Sec. [TT], we discuss the CH- 
FHDG equation and the origin of the concentration-dependent square gradient coefficient 
that has been proposed by KM to cause pinning in off-critical blends. The discretization 
and numerical integration scheme used to solve this equation, and our numerical results, are 
presented in Sec. |T| and discussed in Sec. [TV|. Finally, a summary of our main conclusions, 
and speculations on possible mechanisms of pinning in blends, is discussed in Sec. [V|. 



II. THEORETICAL MODEL 



Model blends are typically described by the Flory-Huggins-De Gennes free energy func- 



tional 19-21 
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where 0(r) is the local concentration of component A (so that 1-0 is the concentration 
of component B), a a and a B are the Kuhn lengths of the two species, A is an effective 



interaction distance between monomers, and the Flory-Huggins free energy is p2| , |23|j 



^ = ± M+ ^ln(l-0 +x «l-0, (2.3) 
k B T N A N B 

where Na {Nb) is the degree of polymerization of chains A (B). Whereas in small molecule 
mixtures the square gradient coefficient is enthalpic (arising from short range interactions 



between molecules) and independent of the local concentration ||24j| , De Gennes proposed 
that the connectivity of polymer molecules in inhomogeneous blends manifests itself through 
an additional, concentration-dependent contribution to the square gradient coefficient k(0) 
fT9|| . The expression for k(0) in Eq. |2]2| was derived to be consistent with the random phase 
approximation result for the inverse structure factor of an incompressible polymer blend 



S 1(q) N A <j> D{q 2 R A ) + N B (1 - 4>o)D(?# B ) 2X ' (2 ' 4) 

Here R 2 is the average square radius of gyration of species i, O is the average value of the 
concentration and the Debye function is D(x) = 2[x — 1 + e~~ x ]/x 2 , with x = q 2 R 2 . In 
the weak segregation limit, the interfacial width is much larger than the chain dimensions 
||24j| , so that the length scales of interest are larger than Ri {q 2 R 2 <C 1), and thus the 
Debye function may be approximated by D~ l (x) = 1 + x/3 + 0(x 2 ). The square gradient 
coefficient in Eq. [Z72] is then obtained from the coefficient of the q 2 term in the Taylor 
expansion of the inverse structure factor, which is related to the free energy functional by 
21| S^q) = (/cbT) -1 ^ 2 ^/^ 2 , where the r.h.s. is evaluated in q-space. 

Because of the approximations made in the calculation of the square gradient expression, 
the Flory-Huggins-de Gennes free energy functional describes the physics of blends in the 
weak segregation limit, and on length scales much larger than the average chain dimension 
fl^ , |2~4] , |2~7|] . In strongly-segregating blends (xN ^> 1) for which x is small but iV — ► oo, 
Eq. 2.2 with a different pref actor in the ^-dependent part is typically used 



The local part of the free energy (Eq. |2.3|) has the same Ginzburg-Landau type of double- 



well structure as small molecules or Ising-like systems . Thus, the only difference between 
the free energy functionals for the simplest small molecule and polymeric systems arises from 
the chain connectivity, which is expressed in the FHDG functional through the reduction 
of the entropic part of the local term, and by the ^-dependent part of the square-gradient 
coefficient. KM proposed that the entropic contribution to the nonlocal part of the free en- 
ergy, namely the concentration-dependent square- gradient coefficient, provides the barrier to 
coarsening which, when combined with the percolation-to-cluster transition, causes pinning 
of the off-critical phase-separating blend. In the next section, we re-examine the numerical 
solution of the time evolution of the CH-FHDG equation for both critical and off-critical 
blends. Specifically, we show that for this model pinning is not observed in the continuum 
limit, regardless of the blend composition, although a dynamical exponent slightly smaller 
than 1/3 is found. 

The theoretical description of spinodal decomposition in binary blends is based on the 
Cahn-Hilliard-Cook equation for the time evolution of the concentration, originally derived 
for small molecule systems ]7|J3T[]: 



|=V.(M W V^j + ,(M). (2.5) 

In this equation, M(0) is the mobility, F{<p} is the coarse-grained free energy functional, 
and 77 is thermal noise. For polymers, the free energy functional is typically taken to be of 
the Flory-Huggins-De Gennes form in Eq. [2.1| , but more general free energy functionals may 
be included. 

In the following, we will always consider for simplicity a symmetric blend, for which 
Na = Nb = iV and o~a = o~b = o~. In this case, the mobility 

MO) = ND(j)(l - 0) (2.6) 

has been proposed ]T9|] , where D is the self-diffusion coefficient. We will also take the effective 
interaction distance A equal to the Kuhn length a. Substituting in Eq. [2.5| the functional 



derivative of F and the expression for M, the time evolution of the concentration is given 
by: 
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180(1-0)/ ' ' 360 2 (1 - 
Note that in Eq. |2.7j we have neglected the thermal noise term. Since we are interested 
in the presence or absence of pinning due to the FHDG free energy functional alone, and 
since it has been shown that the presence of thermal noise in the analogue of this equation 
for small molecule systems does not influence the scaling function or the growth exponent 



during coarsening [|K],[32|], we will neglect noise in our simulations |33|j34 |. 

Eq. pTFI can be rescaled so that it depends on dimensionless space and time variables. 



The transformation, valid only in the unstable region, is the following 
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where Xs = l/(2iV0 o (l — O )) gives the spinodal curve and O is the average value of 
concentration. (This rescaling differs from the rescaling commonly used in experiments by 



a simple numerical factor.) After this transformation, Eq. |2.7| becomes [15 
<90(x,r) 1 



dr 



1 



Xc 



2(X - Xs) 



In 



; X 

'x~Xs 



a 



„ 2 , (1 - 20W 2 .„ ,. 2 
V ^ + 360 2 (1-0) 2 (V ^ 



(2.9) 



III. NUMERICAL SOLUTION 



The solution of this continuum equation ( [2.9| , which was first studied by KM [|T5 
describes the time evolution of the concentration field after a quench to \ > Xs i n the 
unstable region [BBIRII. The initial condition before the quench, corresponding to high 



temperature, is given by a uniform field with random fluctuations about its average value 
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0o- At early times following the quench, the uniform concentration is unstable with respect to 
the long wavelength fluctuations arising from the initial condition, and the two components 
begin to spatially separate. Domains rich in one or the other component form, and then 
coarsen so as to remove interfaces and minimize the free energy. In small-molecule mixtures 
described by, e.g., the Ginzburg-Landau free energy functional, these domains coarsen until 
phase separation is complete regardless of the relative composition and the presence of 
thermal noise. Our goal in this paper is to determine if the same is true for the FHDG free 



energy functional [[34 



In the late stages of decomposition, the system can be characterized by the evolution of 



the typical size of growing domains. The scaling hypothesis p5| states that this length L (as 
calculated from, e.g., the inverse of the peak position of the structure factor, the inverse of 
the first moment of the structure factor, the position of the first zero in the pair correlation 
function, etc.) evolves in time according to 

L ~ r a . (3.1) 

The choice of one particular definition of L is dictated only by convenience. Numerical 
simulations, experiments, and analytical results strongly support the validity of the scaling 
hypothesis in small molecule systems, giving in the absence of hydrodynamic forces the 
value a = 1/3 independent of quench depth, relative composition, and the presence of noise 
While polymers are believed to belong to the same static universality class as small 



molecules, the situation is less clear with respect to dynamics. 

To study the kinetics of spinodal decomposition in polymer blends, we numerically inte- 
grate Eq. |2.9| via a finite difference scheme for both time and space variables. The continuous 
space of position vectors is replaced by n 3 sites on a simple cubic lattice with mesh size (lat- 
tice spacing) Ax. The temporal discretization is achieved by replacing the continuous time 
variable r by a series of m discrete time steps of duration At. The value of the concentra- 
tion field at all sites at each time step is then computed by a first-order Euler numerical 
integration scheme [36|, described in detail in the Appendix. 



Although a large time step and mesh size would speed up the computation, the mesh size 
must be chosen carefully so as to be smaller than the smallest important length scale in the 
problem at all times — here the interfacial width, which decreases in time until the latest 
stages of demixing. The size of the time step is in turn limited by the mesh size. A time 
step that is too large could generate instabilities and spurious solutions |3"2] , |3"6T ] . Thus, these 
discrete variables must be chosen carefully in concert. A linear stability analysis can be 
helpful in suggesting reasonable trial values. If the algorithm is stable with these values of 
Ax and At, one can then vary them to find optimum values and to ensure that the solution 
is accurate and independent of the choice of these parameters. 

We have studied the effect on the numerical solution of Eq. [Z]9] of changing Ax; At 
is changed suitably so as to maintain stability. The boundary conditions are periodic in 
all three directions. Initial conditions are given by random values of the concentration 
field, with average 0o and a flat distribution between O — A0 O and 0o + A0 O - A0 O has a 
strong influence on the behavior during the initial regime, but does not affect the late stage 
behavior. Therefore we fix A0 O — 0.1. Several realizations of the initial conditions were 
averaged together for every set of parameters. 

The phase separation is monitored visually in real space, and quantitatively by deter- 
mining the time evolution of fci(r), the first moment of the spherically-averaged structure 



factor [13]; this is the inverse length that is used to determine the exponent a. For the 



purposes of comparison with previous studies fl5|]3~4]] , we take \ to be related to T (K) by 
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X = 0.326/T-2.3 • 10" 4 , (3.2) 

and fix T c = 62 °C. The system was quenched to temperatures T = 54.5,49 and 25 °C 
for critical composition (0 O = 0.5), and to temperatures T = 35 and 15 °C for off-critical 
composition (<p = 0.4). 

We first consider the solution of Eq. |2.9| obtained with a mesh size Ax = 1, time step 
At = 0.01 and n = 32 (Fig. |l|). This choice of mesh size gives, e.g. for T = 35 °C and 
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0o — 0.4, an equivalent dimensional mesh size of Ar m 6R g , where R g is the average chain 
radius of gyration. With these choices, we are exactly repeating the integration of Eq. |2.9| 
previously performed by Kotnis and Muthukumar. Our results reproduce their findings. For 
critical quenches (0o — 0.5), after an initial transient, the system enters the late stage regime, 
where k\ decays in time with a power law. The exponent a appears to be smaller than 1/3, 
the value expected for spinodal decomposition in small molecule systems and for polymer 
blends in the intermediate stages of demixing (ie. without hydrodynamics). For off-critical 
quenches (0o = 0.4), domain growth stops before the phase separation is complete. 

Results change drastically when the mesh size is reduced to Ax = 0.5, for which the time 
step must be reduced to Ar = 0.002 to maintain numerical stability (Fig. 0). (Note that in 
this case we must take n = 64 to keep the system size the same as before — n ■ Ax = 32). 
This choice of mesh size gives, e.g. for T = 35 °C and 0o — 0.4, an equivalent dimensional 
mesh size of Ar pa 3R g . In this case, the late stage behavior is the same for both critical and 
off-critical quenches, even after the blend undergoes the percolation-to-cluster transition. 
After an initial transient, the late stage scaling regime is reached: k\ decays as a power 
law and no pinning is observed. The values of a can be computed from the slopes of the 
curves in Fig. |2] and are reported in Table |. In all cases a is found to be greater than 1/4 
and smaller than 1/3. We believe that the latter is the true asymptotic value for this free 
energy functional; the small systematic error may be attributed to the crossover from the 
preasymptotic regime and, possibly, to the numerical slowing down given by a still oversized 
mesh size (see Sec. fV\j . Increased accuracy might be obtained by using an even smaller 
value for Ax, but this implies a suitable reduction in the value of Ar, and as a result, a 
prohibitively large computation time. 

The critical relevance of the mesh size to the late stage behavior of the system is also 
evident from Fig. |3], where a "pinned" state obtained with Ax = 1 is "depinned" by reducing 
the mesh size to 0.5. 
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IV. DISCUSSION 



The absence of pinning for small mesh sizes in our simulations even for off-critical 
quenches clearly shows that what is observed for Ax = 1 is not a physical effect but only 



an artifact of the discretization. By integrating Eq. 2^ via a discretization scheme, we are 
actually changing the model under consideration: the solution of the discrete model exactly 
reproduces that of the continuum equation only in the limit where At and Ax approach 
zero. For this reason one must always confirm that the numerical results are independent 
of the values of the discretization variables. In particular, it was recently shown that an 
oversized mesh size can cause a non-physical freezing of interfacial motion for systems with 
non-conserved order parameter |3£| . "Spurious pinning" was already noted by Rogers et al. 



32] in a conserved order parameter system, who showed that for the Cahn-Hilliard equation 



with the Ginzburg-Landau free energy functional, a mesh size Ax > 1.7 causes an unphysical 
decrease in the effective growth exponent. 

In general, the solution of the discrete model should reproduce the behavior of the 
continuum equation if Ax is much less than the smallest physical length modeled in the 
system. In the case under investigation, we are studying a polymer blend in the weak 
segregation limit (% > XcX ^ !)• The smallest physical length that must be resolved is 
the interfacial width, which during the late stages of phase separation is of the order of the 
correlation length £ f24 |: 



which is close to unity in rescaled units. Choosing Ax = 1 implies that we are resolving the 
system at a length equal to the correlation length. The spatial derivatives at the interfaces 
are consequently computed inaccurately with this mesh size, thereby producing unphysical 
results. This picture is confirmed by Fig. f| showing that when Ax = 1 the interface appears 
only one mesh size wide. When Ax = 0.5 the interface is smoother, and larger than the 
mesh size; hence no pinning occurs. This effect also explains the low estimates of a for 
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critical quenches when Ax = 1 — the sharpness of the interface unphysically slows down 
the evolution of the solution, even if it is not sufficient to pin it. Possibly also the results 
for Ax = 0.5 are slightly biased by this effect. 

It is possible to understand the origin of this problem in another way by looking directly 
at the equation of motion. Consider for simplicity a one-dimensional small molecule system 
(i.e. the square-gradient coefficient k and the mobility M are constant), where the concen- 
tration profile goes from one bulk value (0i) at site Xi — Ax to the other (0 2 ) at site Xi + Ax, 
through an interface. For the domain size to grow, the interface must move, and thus (f>(xi) 
must change from <pi to 02- The driving force for this change is the square gradient term in 
the free energy, which yields a Laplacian in the functional derivative of F. This force must 
overcome the double well potential given by the local term in the free energy expression, as 
stated by the Cahn-Hilliard equation: 

| = MV*(g- K VV). (4.2) 

In the discrete version of this equation, the local part does not depend on Ax, while the 
Laplacian is given by: 

V ^ = tw-^ + A;r ) + <K X * - Ax ) - 2< ^)] • ( 4 - 3 ) 

(Ax) z 

When we increase Ax the denominator grows indefinitely, while the numerator is bounded 
above by 02 + 01 — 20(xj) = const. Then, increasing Ax decreases the value of the Laplacian 
term, while the local term is unchanged. For Ax large enough the Laplacian cannot overcome 
the local term, the solution stops evolving, and the system artificially "pins" . 



V. CONCLUSION 

In summary, we have shown that both critical and off-critical polymer blends described 
by the Flory-Huggins-De Gennes free energy functional undergo spinodal decomposition 
via the Cahn-Hilliard equation without pinning of the domain growth as observed in some 
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experiments. Even in the absence of thermal noise |34|, the solution of the discretized 
equation of motion shows coarsening without evidence of pinning, regardless of the relative 
concentration of the blend components. We have also shown that previous solutions of the 
CH-FHDG equation that exhibited pinning were artifacts of an oversized mesh size used in 
the discretization and numerical integration of the equation of motion, and not a result of 



the concentration-dependence of the square gradient coefficient as previously suggested |15 
This suggests the FHDG free energy functional alone, as written in Eq. ^T] is not sufficient 
to describe the physics responsible for pinning in real blends. 

Evidently, a model able to describe the arrested growth observed in experiments must 
include additional physical ingredients. One should consider that experimental blends which 
exhibit pinning often contain components which are not simple homopolymers; such is the 
case with Hashimoto's random copolymer/homopolymer blends, as well as with Hasegawa's 
liquid-crystalline polymer/homopolymer blend, in which the liquid-crystalline component is 
anistropic at the quench temperatures at which pinning was observed. These blends may 
not be describable by the simple Flory-Huggins-De Gennes expression, and consequently 
we should not expect the CH-FHDG equation to mimic their behavior during spinodal 
decomposition. It is also important to note that a significant fraction of homopolymer blends 
never exhibit pinning, regardless of the relative composition. However, for those that do, it 
is possible that either the mobility, or free energy functional, or both, must be appropriately 
modified. The simplified expression for mobility (Eq. [2.6|) used in the CH-FHDG equation, 



which was originally derived only for the special case of a perfectly symmetric blend ||T9 
may neglect important contributions to mobility arising from the connectivity within the 
polymer chains. 

With respect to possible free energy modifications, it is well known that interfacial growth 
can be slowed or stopped by decreasing the interfacial tension. This can be achieved by, 
e.g. using surfactants in the case of small molecules f55[l , diblock copolymers in the case of 
polymers pEP ], or impurities in the case of alloys [ffl |. If the experimental blends contain 



even a small number of impurities, or specific interaction regions that act as impurities 
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42|, these impurities could arrest the demixing process and cause pinning. Finally, the 
FHDG free energy functional in Eq. |2.1| describes incompressible blends; real blends are in 
fact compressible. A coupling of concentration fluctuations and density fluctuations may 



be responsible for pinning in blends ||43|| , in which case a reformulation of the free energy 
functional as well as the addition of a second order parameter field is necessary. 

We are extremely grateful to J. Douglas, B. Hammouda, P. Gallagher, C. Han, E. Di- 
Marzio, G. McFadden, A. Coniglio, F. Corberi and J. Warren, and especially to M. Muthuku- 
mar, A. Chakrabarti and G. Brown, for useful discussions. We thank the Center for Com- 
putational Science at Boston University and the University of Maryland for generous use 
of their CM-5. CC would like to thank the Structure and Mechanics Group in the Poly- 
mers Division at NIST, and the NIST Center for Theoretical and Computational Materials 
Science, for their hospitality. 

VI. APPENDIX 

The numerical solution of Eq. |2.9| is performed via iteration of the following map: 

aim 



This map, given the value of the concentration field 4>™jk a ^ time mAr at each of the 
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sites of a simple cubic lattice with mesh size Ax, yields the value of 4>^k a ^ eacn s ^ e a ^ 
time (m + l)Ar. Note that the mesh size is taken to be the same in all directions. 

The time derivative <90™ k j 'dr is given by the discretization of the left hand side of 
Eq. |2l], with spatial derivatives centrally discretized. This means the chemical potential 
IJ^j k th & t appears in square brackets in Eq. |2]9| is computed using: 

and 

V 0(X ' T) (Axf + {Axf + 
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$j,fc+i 2 ^,fc + <t>T,j,k-i 

(Ax) 2 ' 

The divergence of the product 0(1 — 0)V/x is evaluated as follows: 

/ 7 m _ 7 m \ 



2Ax 



where 



/ ..m ..ui \ 



and Fj™ fc and Z t ™ fc are defined accordingly. 
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FIGURES 

FIG. 1. Plot of log(/ci) vs log(r) for a system with n = 32 and Ax = 1. Each curve is the 
average over 5 realizations with different initial conditions. 

FIG. 2. Plot of log(/ci) vs log(r) for a system with n = 64 and Ax = 0.5. Each curve for 
critical quenches is the average over 5 realizations with different initial conditions. Each curve for 
off-critical quenches is averaged over 10 realizations. 

FIG. 3. (a) Snapshot of a 2D slice of a system of size n = 32 and Ax = 1 at r = 100. (b) The 
same system for r = 200. (c) The same system in (a) at r = 200 when Ax is switched to 0.5 at 
r = 100. The decrease of mesh size allows the system to phase-separate without pinning. 

FIG. 4. Interface profile for (a) Ax = 1 and (b) Ax = 0.5, for n = 32, O = 0.4, T = 15 and 
r = 100. 
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TABLES 



Critical quenches (A(f> = 0.5) 
T 25 49 54.5 

a 0.263 ±0.001 0.286 ± 0.002 0.270 ± 0.001 

Off-critical quenches (A0 O = 0.4) 
T 15 35 

a 0.290 ±0.001 0.306 ± 0.001 

TABLE I. Values of the dynamic exponent a for n = 64 and Ax = 0.5, computed for r > 10. 
The exponents and errors were computed using a linear regression fit of the average values plotted 
in Fig. 2. 
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